function printout_D_rho_of_nonsingular_modes(m, CONSTS, plot_data)

    a     = CONSTS.a;
    eps = CONSTS.eps;
    g = CONSTS.g;
    d = CONSTS.d;
    k0 = CONSTS.k0;

    rho_high = 50;
    num_rho_points = 20000; %200
    rho_step = rho_high/num_rho_points;
    rho_1_vec = (0:rho_step:a)';
    rho_2_vec = (a:rho_step:rho_high)';
    r1_vec = rho_1_vec./a;
    r2_vec = rho_2_vec./a;

    [E_rho_inside,   E_phi_inside,   E_z_inside]  = func_E_inside(r1_vec, m, CONSTS);
    [E_rho_outside,  E_phi_outside,  E_z_outside] = func_E_outside(r2_vec, m, CONSTS);
    
    D_rho_outside = E_rho_outside;
    D_rho_inside = eps*E_rho_inside - 1i*g*E_phi_inside;

    if(plot_data)
        
        figure; plot(r1_vec, real(D_rho_inside),  'b.-', ...
                     r1_vec, imag(D_rho_inside),     'r.-', ...
                     r2_vec, real(D_rho_outside), 'b.-', ...
                     r2_vec, imag(D_rho_outside),    'r.-');
        title('D_{\rho}'); xlabel('\rho/a'); ylabel('D_{\rho}');
        legend('Re D_{\rho}', 'Im D_{\rho}', 'location', 'SouthEast');
        
    end
    
end